function [theta,daligned]=principal_orientations_ext(des,thresh1,thresh2,N,L)
% function principal_orientations_ext(a,N,L)
%
%  calculate principal orientations of a flow field by maximizing its sum
%  of projection coefficients on the basis functions \phi_{k,1}. See the
%  paper for more information.
%
%   INPUTS:
%           des - matrix of descriptors, each COLUMN is a descriptor.
%           thresh2- threshold for pruning the principal orientations. 
%           N -     positive order of the Laurent series
%           L -     negative order of the Laurent series
%                   numel(des)=(N+L+1)*2
%
%
%   See Appendix of our journal paper for more details
%
%
%
% This file is part of the source-code demonstrating the method published 
% in the paper: 
%
% Wei Liu and Eraldo Ribeiro. Scale and Rotation Invariant Detection 
%      of Singular Patterns in Vector Flow Fields. IAPR International 
%      Workshop on Structural Syntactic Pattern Recognition (S-SSPR), 
%      Cesne, Turkey, 2010.
%
% BIBTEX ENTRY:
%
% @inproceedings{WeiRibeiroSSPR2010,
%   author    = {Wei Liu and Eraldo Ribeiro},
%   title     = {Scale and Rotation Invariant Detection of Singular 
%                Patterns in Vector Flow Fields},
%   booktitle = {IAPR International Workshop on Structural Syntactic 
%                Pattern Recognition (S-SSPR)},
%   year      = {2010},
%   pages     = {XXX-XXX}
% }
%
% The source-code can be obtained from 
%        http://www.cs.fit.edu/~eribeiro/flowdetector/
% 
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you 
% maintain the name of the original authors.
%
% (c) 2010 by Wei Liu and Eraldo Ribeiro 
% eribeiro@cs.fit.edu
% Florida Instititute of Technology, Melbourne, Florida, U.S.A.
%
%
%   The key idea is to turn the trigonometric equation into a polynomial
%   one using cos k\theta=(t^k+t^{-k})/2
%
if(nargin<4)
    error('must have 4 inputs');
end

if(numel(des)~=(N+L+1)*2)
    error('dimension of descriptor does not match');
end

kk=[-L-1:N-1]';                                   %k-1;   k=-L, -L+1,...., N
M=max(N-1,L+1);

cc=zeros(2*M+1,1);% the coefficient array, see MATLAB ROOTS function for more details

id1=M+kk+1;id2=M-kk+1;
ak1=des(1:2:end).*kk;
ak2=des(2:2:end).*kk;
cc(id1)=cc(id1)+ak1+ak2*1i;
cc(id2)=cc(id2)+ak2*1i-ak1;

%%
cc=flipdim(cc,1); %flip the coefficient matrix so that the order of the coefficients is descending
t=roots(cc);% solve the polynomial equation
magt=t.*conj(t);
I=abs(magt-1)<1e-4;     %throw away the roots that does not satisfy the trigonometry relationship
t=t(I);
%%
if(numel(t)<1)
    theta=[];
    dsig=[];
    return;
end

%%

%%
sinalpha=real((t-1./t)./2i);
cosalpha=real((t+1./t)./2);
theta=atan2(sinalpha,cosalpha);
%% calculate the second-order derivatives to find the local maxima, and remove
%% the local minima.
kk=kk';
alpha=repmat(theta,[1,numel(kk)]);
tkk=repmat(kk,[numel(theta),1]);
cos_kk_alpha=cos(alpha.*tkk);
sin_kk_alpha=sin(alpha.*tkk);
tdes=repmat(des',[numel(theta),1]);
j2=-tkk.^2.*(tdes(:,1:2:end).*cos_kk_alpha-tdes(:,2:2:end).*sin_kk_alpha);
j2=sum(j2,2);
theta=theta(j2<0);%remove the local minima

%% now calculate \sum a'_{k,1} (\theta) and see if the alignment is
%% significant
cos_kk_alpha=cos_kk_alpha(j2<0,:);
sin_kk_alpha=sin_kk_alpha(j2<0,:);
tdes=repmat(des',[numel(theta),1]);
aligned=sum(cos_kk_alpha.*tdes(:,1:2:end)-sin_kk_alpha.*tdes(:,2:2:end),2);
maligned=max(aligned);
idx=and(aligned>maligned*thresh2,aligned>thresh1);
theta=theta(idx);
daligned=aligned(idx);
return;
